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SOME  WINDOWS  WITH  VERY  GOOD 
SIDELOBE  BEHAVIOR; 
APPLICATION  TO  DISCRETE 
HILBERT  TRANSFORM 


INTRODUCTION 

The  use  of  temporal  weightings  for  spectral  analysis,  with  good  sidelobe  behavior 
and  small  bias,  is  well  established  and  documented  in  Harris  (ref.  1).  However, 
some  of  the  plots  of  the  spectral  windows  are  not  correct  and  do  not  have  the  op¬ 
timal  sidelobe  levels  claimed.  We  will  present  the  corrected  plots  and  some  ad¬ 
ditional  windows  with  optimal  properties.  Finally,  some  of  these  results  will  be 
applied  to  obtain  a  discrete  Hilbert  transform. 

The  temporal  weightings  considered  will  be  continuous  functions  of  time  (except 
possibly  at  t  =  ±  L/2)  and  duration  L;  i.e., 


w(t)  =  0  for  1 1| >L/2  .  (1) 


The  Fourier  transform  of  the  temporal  weighting  is  the  spectral  window 

1/2 

W(f)=  f  dt  w(t)  exp(-i2nft)  (2) 

-1-/2 

and  is  a  continuous  function  of  frequency,  defined  for  all  f.  Notice  the  notational 
convention  adopted  here:  a  weighting  is  applied  multiplicatively  in  one  domain,  and 
its  Fourier-transform  (called  a  window)  occurs  as  a  convolution  in  the  other 
domain. 

All  the  window  results  presented  here  are  obtained  by  exact  analytical  evaluation 
of  (2)  and  are  valid  for  all  values  of  f.  However,  window  (2)  can  be  approximately 
evaluated  at  any  f,  by  means  of  some  numerical  integration  rule  (such  as 
Trapezoidal),  by  choosing  increment  A  =  L/M,  where  M  is  a  large  integer.  These 
latter  results  are  not  adequate  for  |  f  |>.5/A,  because  the  approximation  yielded  by 
this  numerical  integration  procedure  has  period  1/A  in  f.  Furthermore,  if  we  limit 
the  frequencies  f,  at  which  this  numerical  evaluation  is  conducted,  to  the  values 
n/(NA)  (for  n  and  N  integer),  the  results  can  be  realized  as  an  N-point  discrete 
Fourier  transform  (DFT).  Since  the  frequency  spacing  at  which  these  values  occur  is 
(NA)  1  and  the  width  of  the  spectral  window  (2)  is  of  the  order  1/L  =  (MA)  ',  we 
would  also  require  N>M  if  we  desire  to  observe  fairly  closely  the  changes  in  the 
window  (2)  by  means  of  an  N-point  DFT.  There  is  no  fundamental  restriction  on  the 
relative  sizes  of  M  and  N;  however,  M  must  be  large  in  order  to  obtain  an  accurate 
approximation  to  (2). 
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GENERAL  WEIGHTING  CONSIDERATIONS 

The  temporal  weightings  of  interest  here  are  of  the  form 

i  K 

w(t)=  -r-  X  ak  cos(2?rkt/L)  for  |t|<L/2,  (3) 

L  k  =  0 

where  {ak}g  are  real  constants.  The  weighting  is  symmetric  about  t  =  0  and  possesses 
all  orders  of  derivatives  for  |  t  |<L/2;  however,  discontinuities  in  w(t),  defined  by 
(1)  and  (3),  or  in  its  derivatives,  occur  at  t  =  ±  L/2.  These  discontinuities  dictate 
the  asymptotic  behavior  for  large  j  f  |  of  W(f)  in  (2).  Without  loss  of  generality,  the 
weighting  is  normalized  according  to 

K  K 

X  ak=  1;  w(0)=-j-  2  ak  =  4- .  (4) 

k  =  0  L  k  =  0  L 

Observe  from  (3)  that 

i  1  K 

w(±^)  =  lim  w(t)=-p  X  (-l)k  ak  ,  (5) 

z  |l|— L/2-  L  k  -0 

which  may  or  may  not  be  equal  to  zero.  If  (5)  is  not  zero,  then  weighting  w(t)  is 

discontinuous  at  t  —  ±  L/2,  and  window  W(f)  will  decay  only  as  I/f  for  large  [  f  f. 

However,  if  (5)  is  zero,  then  w(t)  is  continuous  for  all  t.  Also,  w  '(t)  is  continuous 
for  all  t,  since  we  always  have 

,  K 

w'(t)=  -  T  kak  sin(2nkt/L)  for  jt|<L/2  ,  (6) 

L  k  =  0 

and 


lim  w'(t)  =  0;  w'(t)  =  0  for  |t|>L/2  .  (7) 


The  last  property  follows  from  (1).  Thus,  when  weighting  values  w(  ±  L/2)  in  (5)  are 
zero,  w(t)  and  w  '(t)  are  both  continuous  for  all  t. 

However,  w  "(t)  may  then  not  be  continuous  at  t  =  ±  L/2.  We  have  from  (6), 


and 


w"(t)  = 


4n2 

LJ 


K 


X  k2ak 

k  --  0 


cos(2nkt/L) 


for 


1 1 1  <L/2 


A  2  ^ 

lim  w"(t)  =  -  X  (-l)kk2ak, 
UI-I./2-  L  k=0 


(8) 


(9) 
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which  may  or  may  not  be  zero.  If  (9)  is  not  zero,  then  w”(t)  is  discontinuous  at 
t  =  ±  L/2,  and  W(f)  will  decay  as  1/f 3  for  large  |  f  |. 

Conversely,  if  (9)  is  zero,  then  w''(t)  is  continuous  for  all  t,  and  it  follows 
(similarly  to  above)  that  w"'(t)  is  continuous  for  all  t.  Then  W(f)  decays  at  least  as 
fast  as  1/f5  for  large  |  f  | .  We  will  have  occasion  to  use  these  relations  later. 

The  spectral  window  corresponding  to  w(t)  in  (1)  and  (3)  is  given  by  (2)  as  the 
closed  form  expression 


W(f)=  —  sin(nLf)X  — —  for  all  f; 

k  -0  LTP-k2 


W(n/L)  = 


a0,  n  =  0 

}a!n|,n#0 

k  \  2 


( 1 0  A ) 


If  we  expand  (L2f 2-k2)“ 1  in  a  power  series  in  (  q)  ",  we  obtain 


W(f)  = 


1  ,r4.  I  (-Dk  k2m  ak  for  |  L f |  > K  . 


sin(nl  Q 
nLf  <n~ o (L02nl 


(10B) 


Thus,  the  quantities  (5)  and  (9)  considered  above  are  simply  the  m  =  0  and  m=  1 
coefficients  of  this  expansion;  the  asymptotic  behavior  of  (10B)  depends  on  the  first 
non-zero  term  in  the  m-series,  and  will  be  plotted  in  the  following  figures  as  a  dotted 
line. 
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DISCRETE  FOURIER  SERIES 


When  continuous  temporal  weighting  is  used  in  digital  processing,  it  is  sampled 
and  often  transformed  into  the  frequency  domain,  where  its  effect  can  be  included 
as  a  convolution  of  its  discrete  Fourier  series  with  the  data  DFT.  To  evaluate  this 
discrete  Fourier  series,  we  begin  by  delaying  the  temporal  weight  to  the  interval 
(0,  L): 


wD(t)sw(t  ~y)  =  X  X  (-l)kakcos(2nkt/L) 

K  k 

=  -j-  X  (-l)k  £k  a,k I  exp(i2nkt/L)  for  (Kt<  L, 

L  k  =  -K 


where 


1 ,  k  =  0 
k  */2 ,  k*0 


(ID 


(12) 


Let  the  sampling  interval  on  weighting  wD(t)  be 

A  =  L/M  ,  (13) 

where  M  is  even;  then  samples  (temporal  weights) 

1  ^ 

wD(mA)  =  -r-  (-l)k£k  aiki  exp(i2rckm/M)  for  (Km<M  .  (14) 

L  k  -K 


Then  for  M>2K,  the  discrete  Fourier  series  is  given  by  the  M-point  DFT  (ref.  2,  ch. 
3)* 


M-l 

£  Aw[:)(mA)  exp(-i2jimn/M)  =  (- 1  )n  en  a,n  for  |n|<M/2  . 

m  -  0 


(15A) 


Thus,  the  effects  of  temporal  weighting  (14)  can  be  incorporated  in  a  digital 
processing  application  as  a  frequency  convolution  of  the  data  DFT  with  the 
sequence 

.  .  .  ,0,0,(-l)K^-aK,  .  .  .,-ya„a0,-yal . (-1)K  yaK,0,0 .  (15B) 

This  is  one  of  the  main  reasons  for  employing  weightings  of  the  form  (14)  in  digital 
processing  applications;  namely,  the  effects  of  temporal  weighting  are  easily  in- 


*  This  consideration  is  different  from  that  mentioned  under  (2), where  we  were 
interested  in  approximately  evaluating  window  W(f)  in  (2)  by  means  of  a  DFT. 
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eluded  by  means  of  convolution  in  the  discrete  frequency  domain  with  a  short 
sequence  of  length  2K  +  1 . 

The  effect  of  sampling  continuous  weighting  w(t)  at  increment  A,  in  so  far  as  the 
effective  window  is  concerned,  is  as  follows.  The  effective  window  is* 


We(f)  =  J  dt  AdA(t)  w(t)  exp(-i2rrft) 

oo 

=  <5,(0  @W(0=  2  W(f- -jj-)  ,  (16) 

4  n  =  -oo 

where  <5a(x)  is  an  infinite  unit-area  impulse  train  in  x  at  spacing  a.  Thus  there  are 
periodic  replications  in  We(f)  at  multiples  of  1/A;  the  aliasing  at  f  =  ±.5/A  is  ob¬ 
vious.  All  the  following  results  correspond  to  A  =  0,  i.e.,  continuous  weightings; 
thus,  there  is  no  aliasing. 


Integrals  without  limits  are  over  the  range  of  non-zero  integrand. 
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HARRIS’ WINDOWS 

The  first  weighting  to  be  considered  is  Hanning,  for  which  there  are  only  two 
non-zero  coefficients  in  (3): 


ao=  -5, 


a,=  .5,  V/(0  = 


sin(nLf) 

27tLf(l-L2f2) 


(17) 


Its  power  response  in  dB  is  plotted  versus  Lf  in  figure  1,  normalized  to  the  peak 
response  at  f  =  0;  that  is 


dB  =  10  log|  W(f)/W(0)|2  (18) 

is  plotted.*  The  largest  sidelobe  is  -31 .47  dB.  Since  (5)  is  zero  for  the  weights  in  (1 7), 
but  (9)  is  not,  the  Hanning  window  has  an  asymptotic  decay  of  18  dB/octave;  this 
decay  is  the  first  non-zero  term  in  (10B). 

The  second  weighting  is  Blackman  (ref.  1 ,  eq.  32): 

a()  =  .42,  a,  -  .50,  a2  =  .08  .  (19) 

The  window  is  depicted  in  figure  2.  Again,  (5)  is  zero,  but  (9)  is  not;  thus  the 
asymptotic  decay  is  1 8  dB/octave.  The  largest  sidelobe  is  -58. 1 1  dB. 

The  “exact”  Blackman  weights  are  (ref.  1 ,  p.  63) 

a0  =  7938/18608,  a,  =  9240/18608,  a,  =  1430/18608  .  (20) 

Now  (5)  is  not  zero;  therefore,  the  window  decays  at  only  6  dB/octave  as  shown  in 
figure  3.  However,  the  largest  sidelobe  is  -68.24  dB,  not  -51  dB  as  cited  in  ref.  1, 
fig.  23.  Also,  the  sidelobes  in  figure  3  are  about  6  dB  lower  than  those  reported  in 
ref.  1 ,  fig.  23. 

The  following  four  windows  are  listed  in  the  table  on  page  65  of  ref.  1.  The 
“minimum”  3-term  weights  are 

a0  =  .42323,  a,  =  .49755,  a2  =  .07922  .  (21) 

Since  (5)  is  not  zero,  the  window  decay  is  only  6  dB/octave.  as  shown  in  figure  4. 
The  maximum  sidelobe  is  -70.83  dB,  not  -67  dB  as  reported  in  ref.  1,  figure  24  and 
page  64. 

*  When  the  weighting  is  applied  instead  in  the  lag  domain,  as  for  Blackman-Tukey 
spectral  analysis,  rather  than  in  the  time  domain  as  presumed  here,  the  window 
appears  linearly  rather  than  as  its  square.  In  (his  case,  the  square  must  be  removed 
from  the  definition  in  (18),  and  all  the  plots  require  that  the  dB  numbers  on  the 
ordinate  be  halved.  For  example,  the  peak  Hanning  sidelobe  is  then  - 15.73  dB. 
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The  3-term  weights  are 

a0  =  .44959,  a,  =  .49364,  a2  =  .05677  .  (22) 

The  corresponding  window  is  given  in  figure  5  and  has  a  maximum  sidelobe  of 

-62.05  dB,  rather  than  the  -61  dB  reported  in  ref.  1,  p.  65.  Since  (5)  is  not  zero,  the 
decay  is  only  at  6  dB/octave,  as  indicated  by  the  dotted  line. 

The  “minimum”  4-term  weights  are* 

a0  =  .35875,  a,  =  .48829,  a,  =  .14128,  a3  =  . 01168.  (23) 

Although  (5)  is  not  zero,  it  is  nearly  so.  Therefore,  the  initial  decay  of  the  window  is 
greater  than  6  dB/octave;  however,  it  must  eventually  decay  only  as  6  dB/octave. 
The  maximum  sidelobe  of  the  window  is  indicated  in  figure  6;  it  is  -92  dB,  as 
reported  in  ref.  1 . 

The  4-term  weights  are 

a0  =  .40217,  a,  =  .49703,  a,  =  .09892,  a,  =  .00188  .  (24) 

The  asymptotic  decay  is  only  6  dB/octave,  as  shown  in  figure  7,  and  the  largest 
sidelobe  is  -74.39  dB,  as  claimed  in  ref,  1 ,  p.  65. 


*  These  are  not  the  weights  actually  listed  in  ref.  I,  page  65.  However,  the  values 
listed  there  do  not  add  up  to  1;  accordingly,  we  have  modified  them  according  to  the 
comment  under  eq.  34,  and  made  them  sum  to  I  (with  the  lowest  sidelobe  possible 
by  modifying  just  one  of  the  last  digits). 
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Figure  3.  Exact  Blackman  Window  for  a0  =  7938/ i 8608, 
a,  =9240/18608,  a2  = 1430/18608 


Figure  4.  “Minimum”  3-Term  Window  for  a0=  .42323, 
a ,  =  . 49755,  a2  =  . 07922 


RAPIDLY  DECAYING  WINDOWS  WITH  MINIMAL  SIDELOBES 


It  was  observed  earlier  that  the  window  W(f)  decays  fairly  rapidly  for  large  f  if  (5) 
is  zero,  and  very  rapidly  if  (9)  is  zero.  Such  windows  will  lead  to  spectral  estimates 
that  are  immune  to  strong  interferences  at  frequencies  removed  from  those  of  in¬ 
terest.  In  this  section,  we  will  consider  this  class  of  windows  in  terms  of  the  peak 
sidelobe  and  asymptotic  decay;  the  mainlobe  width  of  each  window  is  not  discussed, 
but  is  easily  determined  from  the  plots.  Discontinuous  weightings  will  be  taken  up 
later. 

If  only  two  coefficients  in  weighting  (3)  are  non-zero,  satisfaction  of  (4),  and 
setting  (5)  equal  to  zero,  yield 


a0  +  a,  =  1,  a0- a,  =0  .  (25) 

The  Hanning  window  satisfies  these  requirements  and  is  plotted  in  figure  1.  Con¬ 
volution  sequence  (15B)  is  simply  -1/4,  1/2,  -1/4. 

Moving  on  to  three  non-zero  coefficients  in  (3),  if  we  satisfy  (4),  and  set  (5)  and 
(9)  equal  to  zero,  we  find 


a0  =  3/8,  a,  =  4/8,  a2  =  1/8  .  (26) 

The  weighting  is 

w(t)  =  cos4(nt/ L)  for  |t|<L/2  .  (27) 

From  (15B),  the  discrete  Fourier  series  for  convolution  is 

1,-4,  6,  -4,  1  , 

16  (28) 

which  are  simply  the  binomial  coefficients.  As  noted  under  (9),  since  the  third 
derivative  of  w(t)  is  continuous  at  t  =  ±  L/2,  the  window  decays  at  a  30  dB/octave 
rate.  The  plot  in  figure  8  indicates  that  the  largest  sidelobe  is  -46.74  dB. 

Instead  of  forcing  (9)  equal  to  zero,  we  can  use  the  one  degree  of  freedom  left, 
after  (4)  is  satisfied  and  (5)  is  set  equal  to  zero,  to  minimize  the  maximum  sidelobes. 
The  optimal  weights  are  found  to 

a0  =  .40897,  a,  =  .5,  a2  =  .09103  .  (29) 

The  corresponding  window  is  presented  in  figure  9.  The  asymptotic  decay  is  18 
dB/octave,  and  the  two  equal  sidelobes  are  of  size  -64. 19  dB.  This  is  6. 1  dB  better 
than  the  -58.1  dB  sidelobe  of  the  Blackman  window,  yet  the  asymptotic  decays  are 
equal.  Although  the  maximum  sidelobe  of  the  “minimum"  3-term  window  in  figure 
4  is  6.6  dB  better,  that  decay  is  only  6  dB/octave  rather  than  the  18  dB/octave  decay 
here. 


12 


TR  6239 


When  we  consider  four  non-zero  coefficients  in  (3),  we  have  several  alternatives. 
If  we  satisfy  (4),  set  (5)  and  (9)  both  equal  to  zero,  and  also  set  the  fourth  derivative 
of  w(t)  equal  to  zero  at  t  =  ±  L/2,  we  have  four  equations  in  four  unknowns,  with 
solution 


? 

< 


3q,  3| ,  a2,  a2 


10,  15,6,  1 
32 


The  weighting  is 

w(t)=  —  cos6(nt/L)  for  |(|  < 
and  from  (15B),  the  discrete  Fourier  series  for  convolution  is 


-1, 6,-15,  20,  -15,6,  -1 
64 


(30) 


(31) 


(32) 


which  are  again  the  binomial  coefficients.  The  window  decays  at  a  very  fast  rate  of 
42  dB/octave,  since  the  fifth  derivative  of  w(t)  is  continuous  for  all  t.  The  plot  in 
figure  10  shows  the  maximum  sidclobe  to  be  -60.95  dB. 

If  we  satisfy  (4),  and  set  (5)  and  (9)  both  equal  to  zero,  but  use  the  remaining 
degree  of  freedom  to  minimize  the  maximum  sidelobes,  the  optimal  weights  are 
determined  to  be 


a0  =  .338946,  a,  =  .481973,  a2  =  .161054,  a,  =  .018027  .  (33) 

The  window  is  given  in  figure  1 1  and  has  two  equal  sidelobes  of  -82.60  dB.  1  he 

asymptotic  decay  is  39  dB/octave,  since  the  third  derivative  of  w(t)  is  continuous  tor 
all  t.  Comparison  with  the  "minimum”  4-term  window  in  figure  6  reveals  a  dil- 
ference  of  9.4  dB  in  the  maximum  sidelobe;  however,  the  decay  of  figure  1 1  is  much 
faster  at  a  30  dB/octave  rate.  As  far  as  the  4-term  window  in  figure  7  is  concerned, 
figure  II  has  an  8.2  dB  better  maximum  sidelobe  and  a  much  better  decay,  30 
dB/octave  versus  6  dB/octave. 

Finally,  if  we  satisfy  (4),  and  set  only  (5)  equal  to  zero,  and  use  the  remaining  two 
degrees  of  freedom  to  minimize  the  maximum  sidelobes,  the  optimal  weights  are 
found  to  be 

a0  =  .355768,  a,  =  .487396,  a2  =  .144232,  a,  =  .012604  .  (34) 

The  window  is  show  n  in  figure  12  and  has  three  equal  sidelobes  of  -93.32  dB.  Notice 

that  this  level  is  better  than  the  purported  "minimum"  4-term  level  of  -92  dB 
claimed  in  ref.  I,  pp.  64-65;  and  the  asymptotic  decay  is  18  dB  octave,  not  6 
dB/octave.  Furthermore,  this  level  was  achieved  under  the  constraint  of  setting  (5) 
equal  to  zero.  If  we  were  to  eliminate  this  constraint  of  a  continuous  weighting 
function,  a  sidelobe  level  lower  than  -9.7.32  dB  can  be  achieved.  (This  problem  and 
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the  determination  of  the  true  minimum  3-term  window  (to  replace  figure  4)  are 
undertaken  in  the  next  section.)  Comparison  with  the  4-term  window  of  figure  7 
reveals  an  18.9  dB  peak-sidelobe  improvement  in  figure  12  and  a  better  decay  of  18 
dB/ octave  instead  of  6  dB/octave. 


Figure  8.  Window  for  a0=  .375,  a(  =  .5,  a2  =  .125 
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Figure  9.  Window  for  a0  =  .40897,  a]  =  .5,  a2  =  .09103 


Figure  10.  Window  for  ag  =  10/32,  aj  =  15/32,  a2  =  6/32  a3  =  1/32 
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Figure  11.  Window  fora0=  .338946,  a,  =  .481973,  a2  =  .161054,  a3  =  .018027 


Figure  12.  Window  for  a0=  .355768,  a,  =  .487396,  a2=  .144232,  a3=  .012604 


MINIMUM  SIDELOBE  WINDOWS 


If  only  two  coefficients  in  weighting  (3)  are  non-zero,  and  we  disregard  the 
continuity  requirement  (5),  the  one  degree  of  freedom  left,  after  normalization  (4)  is 
satisfied,  can  be  used  to  minimize  the  maximum  sidelobes.  The  result  is  the  familiar 
Hamming  window,  plotted  in  figure  13,  with  coefficients 

a0  =  .53836,  a,  =  .46164  .  (35) 

The  two  equal  peak-sidelobes  are  -43.19  dB,  and  the  asymptotic  decay  is  only  6 
dB/octave,  as  dictated  by  ( 10B)  when  (5)  is  not  zero. 

For  three  non-zero  coefficients  in  (3),  satisfaction  of  (4)  leaves  two  degrees  of 
freedom.  These  can  be  used  to  realize  the  minimum  3-term  window  in  figure  14,  for 
which  the  optima!  coefficients  are 

a0  =  .4243801,  a,  =  .4973406,  a2  =  .0782793  .  (36) 

There  are  three  equal  peak-sidelobes  of  -71.48  dB,  which  is  .65  dB  better  than  figure 
4,  with  an  asymptotic  decay  of  6  dB/octave  for  both. 

When  four  coefficients  are  non-zero  in  (3),  there  are  three  degrees  of  freedom  left 
after  normalization  (4).  The  minimum  4-term  w  indow  results  for  coefficients 

a(l  =  .3635819,  a,  =  .4891775,  a:  =  .1365995,  a,  =  .0106411  ,  (37) 

and  is  shown  in  figure  15.  The  four  equal  sidelobes  are  at  level  -98.17  dB,  which  is 
6. 16  dB  better  than  figure  6,  with  an  asymptotic  decay  of  6  dB/octave  for  both.  This 
sidelobe  level  is  4.85  dB  better  than  figure  12,  but  the  decay  in  figure  12  is  18 
dB/octave. 
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Figure  14.  Minimum  3-Term  Window  for  a0  =  .4243801. 
a ,  =  .4973406.  a,  =  .0782793 
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COMPARISON  WITH  KAISER-BKSSKI,  ANI)  VAN  DKR  MAAS  WINDOWS 

Hie  windows  in  figures  13,  14,  and  15  are  very  similar  to  i he  Kaiser-Bessel 
window.  Specifically,  the  Kaiser-Bessel  weighting  and  window  arc 


w(t)  =-J-  l„  (B\/l  -(2t/L)2  )  for  |t|  <  L/2 


sin  (\/  n-L-f2  -  B:  ) 


for  all  f. 


respectively,  where  B  is  a  parameter.  If  we  choose  B  to  make  the  first  null  of  the 
Kaiser-Bessel  window  lie  at  the  three  alternatives  of  Lf  =  2,3,4  (as  in  figures  13,  14, 
and  15,  respectively),  we  obtain  the  plots  in  figures  16,  17,  and  18.  The 
corresponding  mainlobe  shapes  are  indistinguishable,  and  the  asymptotic  decays  are 
all  6  dB/ octave.  The  immediate  sidelobes  of  the  Kaiser-Bessel  windows  are  several 
dB  larger  than  the  minimum  results  in  figures  13,  14,  and  15,  but  the  distant 
sidelobes  of  the  Kaiser-Bessel  windows  are  over  10  dB  lower  for  the  examples 
considered.  Thus,  a  trade-off  exists  between  the  peak  sidelobe  and  the  distant 
sidelobe  level. 

The  windows  here  are  also  similar  to  the  ideal  impulsive  van  der  Maas  window 
given  by 


w(t)  = 


_B  1|  (Bn/1  -  (2t/L)2  ,JL 


'\  -  (2t/L)2 


+  2  ^  0  “"j")  +  2  ^  +  ~2  ^  *or '' 1 1  ^~2  ' 


W(f)  =  cos  (\/(  nLf )’  -  B-  )  all  f  •  (39) 

This  window  is  characterized  by  having  the  narrowest  possible  mainlobe  width  for  a 
specified  sidelobe  level,  and  vice  versa.  However,  the  window  does  not  decay  at  all 
for  large  f.  The  peak  to  sidelobe  voltage  level  is  SL=cosh  (B),  and  the  first  null  of 
the  window  occurs  at 


-•  (!•") 


Thus,  the  first  null  location  can  be  expressed  in  terms  of  the  sidelobe  level  SL  ac¬ 
cording  to 


<»=(h  [MsLi^un yy 
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Table  1  shows  this  null  location  and  the  actual  location  for  several  of  the  windows 
presented  earlier,  when  the  peak  sidelobes  are  equal;  the  agreement  is  close, 
especially  for  those  windows  with  a  6  dB/octave  decay,  figures  13,  14,  and  15. 


Table  1.  Comparison  of  Null  Location  with  van  der  Maas  Case 


Figure  Number 

Null  Location 

van  der  Maas 

Null  Location 

9 

3 

2.62 

1 1 

4 

3.29 

12 

4 

3.67 

13 

2 

1.87 

14 

3 

2.88 

15 

4 

3.85 

Figure  16.  Kaiser-Bessel  Window  with  First  Null  at  I.f  =  2 
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DISCRETE  HILBERT  TRANSFORM 

For  a  continuous-time  real  waveform  x(t),  the  (voltage-density)  spectrum 


X(f)  =  J  dt  x(t)exp(-i2irf() 


(42) 


is  conjugate  symmetric  about  f  =  0;  i.e.,  X(-f)  =  X*(f).  If  we  suppress  the  negative- 
frequency  components  and  double  the  positive-frequency  components,  we  obtain 
the  analytic  waveform: 


X  ,  (f)  =  [  1  +  sgn(f)  ]  X(f)  , 


and 


x  .  (t)  =  x(t)  +  ixH(t)  . 


Here  xH(t)  is  the  H i  1  bert  transform  of  x(t): 

xH(t)  =  sr®x(D  =  ljdT7~ 


(43) 


(44) 


(45) 


When  we  consider  instead  a  discrete-time  real  waveform  with  samples  occurring 
every  A  seconds,  the  spectrum  has  period  )/A,  and  the  discrete  Hilbert  transform 
(DHT)  takes  another  form  (ref.  2,  ch.  7).  Specifically,  the  Fourier  coefficients  of  an 
N-point  DFT  of  the  waveform  samples  are  multiplied  by  the  sequence  {un|^  ’, 
w  here  (ref.  2,  p.  354) 


u 


1,  n  =  Oand  n  =  N/2 

2,  0  <  n  <  N/2 

0,  N/2  <  n  <  N- ! 


(46) 


However,  this  ideal  characteristic  exists  only  in  a  theoretical  sense;  its  effect  must  be 
approximated  in  practice  (ref.  2,  p.  36)). 

A  method  of  obtaining  an  accurate  DHT  for  long  data  sequences  is  now 
presented.  We  begin  by  considering  a  continuous  cosine  wave  al  frequency  f„: 

x(t)  =  Acos(2rrf0t  +  9)  .  (47) 

The  spectrum  of  this  waveform,  from  the  infinite  transform  (42),  is  a  pair  of  im¬ 
pulses: 


X(f)  =  y  exp  (i0)  d(f-f„)  +  y  exp(-i0)d(f+  f<()  . 


(48) 
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Then,  as  above, 


X,(f)  =  Aexp(ifl)d(f-g  , 
x  +  (t)  =  A  exp  (i2nf0t  +  i 9) 

=  A  cos(2nf0t  +  6)  +  i  A  sin(2nf()t  +  6)  .  (49) 

Thus,  the  Hilbert  transform  of  x(t)  in  (47)  is  the  sine  component. 

The  procedure  above  succeeded  because  we  completely  eliminated  the  negative- 
frequency  components  of  (48).  However,  suppose  we  compute  the  spectrum  from  a 
weighted  time-limited  segment  of  x(t)  in  (47).  Thus: 


y(t)  =  w(t)x(t)  , 

Y(f)  =  W(f)  ®  X(f)  . 

=  y  exp(i0)  W(f-f„)  +  y  cxp(-i0)  W(f+ f„)  ,  (50) 

w  here  weighting  w(t)  satisfies  (1).  Figure  19  illustrates  spillover  of  the  negative-  and 
positive-frequency  components  across  the  f  =  0  border.  This  spillover  can  be 
controlled  by  choosing  a  window  W(f)  with  low  near-by  sidelobes  and  rapid  decay 
of  distant  sidelobes.  Then,  if  Lt'(l  is  large  enough,  when  we  multiply  Y(f)  by 
1  +  sgn(f)  to  yield  Z(f),  there  will  be  little  distortion  from  the  desired  function 

D(f)  =  A  exp(i0)  W(f-f0)  for  all  f  ; 

d(t)  =  Acxp(i2nf0t  +  i0)w(t)  .  (51) 

We  can  then  inverse  Fourier  transform  Z(f)  to  yield  z(t),  and  then  divide  z(t)  by 
weighting  w(t),  in  an  attempt  to  approximate  analytic  waveform  x  4  (t).  This  method 
will  succeed  in  the  interior  region  of  the  particular  segment  of  the  real  data 
waveform  selected,  but  will  be  in  error  near  the  edges  of  the  segment  because  the 
neighboring  data  points  (outside  the  selected  segment)  were  never  accounted  for. 
This  insufficiency  can  be  alleviated  by  using  tapered  temporal  weighting  and  (say) 
50%  overlapped  segments,  but  retaining  only  the  central  50%  of  the  output  data 
points  of  each  segment,  after  the  processing  indicated  above  is  carried  out.  How 
successful  the  technique  is  depends  critically  on  Lf0  and  the  window  W(f). 

For  a  discrete-time  waveform,  the  considerations  arc  similar  to  those  presented 
above,  e  .ccpt  that  we  now  employ  the  discrete  multiplicative  window  (46)  (instead 
of  1  +  sgn(f)),  and  that  Nyquist  frequency  (2A)  1  must  also  be  considered,  in  ad¬ 
dition  to  zero  frequency,  in  terms  of  spillover.  The  procedure  to  be  investigated  is  as 
follows:  a  long  real  data  sequence  Jxn}  is  broken  up  into  segments  of  length  N, 
where  adjacent  segments  overlap  by  50%.  Temporal  weighting  {wn}^  'is 
multiplicatively  applied  to  the  N  data  points  of  a  segment,  and  an  N-point  OFT  is 
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taken  of  this  product  sequence.  Sequence  {u,,}*  1  in  (46)  is  then  applied 
multiplicatively  to  the  complex  DFT  output.  An  inverse  DFT  of  this  product  is 
evaluated,  and  the  result  divided  by  weighting  { wn}^ '.  Then  only  the  interior  50% 
of  these  N  complex  values  are  retained,  and  the  other  50%  discarded.  By  then 
abutting  values  of  adjacent  segment  outputs,  we  obtain  the  analytic  waveform 
sequence  corresponding  to  the  original  real  data  sequence  {xn}.  Very  long  data 
sequences  can  be  handled  in  this  manner;  only  the  extreme  edges  are  subject  to 
significant  error. 

If  we  look  at  some  of  the  better  windows  presented  earlier,  for  example  figures 
11,  12,  and  15,  we  see  that  very  deep  sidelobes  are  realized  and  maintained  for 
|Lf|>4.  Thus  for  a  sinusoid  at  frequency  f0,  spillover  is  small  in  figure  19  if 


4  1  4 

L^'°^2A  L 


(52) 


the  upper  bound  is  due  to  spillover  at  the  Nyquist  frequency.  Since  segment  length  L 
is  composed  of  N  data  points,  (52)  can  be  expressed  as 


4  1  4 

N<f^<2  -  N 


(53) 


For  example,  a  1024-point  DFT  yields 


,004<f()A<  .496  .  (54) 

Thus,  the  majority  of  frequency  components  in  sequence  {xn|  are  capable  of  very 
accurate  Hilbert  transformation.  To  widen  the  bounds  in  (52),  we  must  increase  L 
and  decrease  A.  If  A  is  chosen  to  minimize  aliasing  in  the  first  place,  the  only 
recourse  is  to  increase  N,  the  number  of  points  in  each  data  segment.  Thus,  larger- 
si/.e  DFTs  should  perform  better  in  terms  of  lower  error  of  the  DHT  and  lower 
allowed  frequencies  in  the  given  sequence  {xnJ. 
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SIMULATION  RESULTS 


The  fifteen  weightings  given  earlier  are  summarized  in  table  2.  They  have  been 
tried  as  candidates  for  the  DHT  technique  presented  in  the  previous  section.  Input 
signals  composed  of  real  unit-amplitude  pure  tones  at  several  different  frequencies 
and  phases  have  been  used  to  test  the  accuracy  of  the  DHT  technique.  In  particular, 
tone  frequencies  f0,  taking  the  four  values 

f0A  =  {V^/5,  V^/27,  V2/65,  \6/350}  =  { .283,  .052,  .022,  .004}  (55) 

were  selected.  (The  square  root  was  used  to  avoid  taking  synchronized  periodic 
samples  of  the  input  waveform.  It  was  found  that  samples  synchronized  with  the 
tone  frequency,  such  as  f0A=.25,  gave  inordinately  small  errors.)  The  phases  of 
each  tone  were  taken  as  0  and  -n/2  corresponding  to  a  cosine  and  sine  waveform, 
respectively.  Intermediate-phase  tones  are  available  as  linear  combinations  of  these 
two  cases. 

Table  2.  Window  Characteristics 


Weighting 

Peak  Sidelobe 
<dB) 

As)mplolic  l)ecav 
(dH/ octave) 

I. 

Hanning 

-31.47 

18 

2. 

Blackman 

-58.11 

18 

3. 

Exact  Blackman 

-68.24 

6 

4. 

“Minimum”  3-Term 

-70.83 

6 

5. 

3-Term 

-62.05 

6 

6. 

“Minimum”  4-Term 

92.01 

6 

7. 

4-Term 

-74.39 

6 

8. 

3-Term  with  Continuous  Third  Derivative 

-46.74 

30 

9. 

3-Term  with  Continuous  First  Derivative 

-64.19 

18 

10. 

4-Term  with  Continuous  Fifth  Derivative 

-60.95 

42 

11. 

4-Term  with  Continuous  Third  Derivative 

-82.60 

30 

12. 

4-Term  with  Continuous  First  Derivative 

-93.32 

18 

13. 

Hamming 

-43.19 

6 

14. 

Minimum  3-Term 

-71.48 

6 

15. 

Minimum  4-Term 

-98.17 

6 

Since  the  Hilbert  transform  of  a  cosine  with  arbitrary  phase  is  a  sine  of  the  same 
phase  (see  (49)),  we  can  compute  the  error  of  the  DHT  technique  very  easily  for  pure 
tones.  For  the  retained  interior  50%  output  data  points  of  each  segment,  four  error 
measures  were  evaluated.  They  were  the  maximum  magnitude-error  and  the  average 
magnitude-error  for  both  the  real  and  imaginary  parts  of  the  analytic  waveform 
sequence.  The  two  error  measures  (maximum  and  average)  for  the  real  part  were 
virtually  zero,  because  of  the  trigonometric  function  error  and  the  round-off  noise 
of  the  12-digit  calculator  employed.  The  average  magnitude-error  measure  for  the 
imaginary  part  of  the  analytic  sequence,  i.c.,  the  error  of  the  discrete  Hilbert 
transform,  is  presented  in  tables  3-6  for  the  four  tone  frequencies  listed  in  (55).  A 
program  for  the  DHT  technique  and  the  error  calculation  is  presented  in  the  ap¬ 
pendix. 


The  first  fact  to  observe  from  tables  3-6  is  that  the  quality  of  the  DHT  varies 
greatly  with  the  window.  In  fact,  the  ratio  of  errors  can  become  as  large  as  300,000 
when  we  consider  different  windows  for  a  specified  waveform.  In  order  to  present 
the  numerical  values  more  conveniently  for  comparison,  the  errors  have  been  scaled 
by  109  in  tables  3-5  and  by  106  in  tabic  6.  Thus,  some  of  the  errors  listed  in  tables  3-6 
are  small  indeed. 

The  maximum  magnitude-error,  for  all  the  input  frequencies,  phases,  and 
windows  considered,  was  approximately  4-to-6  times  larger  than  the  average 
magnitude-error.  (It  has  not  been  tabulated  here,  in  order  to  save  space.) 

When  the  data  frequency,  f0,  is  near  the  center  of  the  zero-to-Nyquist  frequency 
band,  the  best  windows  are  observed  from  tables  2  and  3  to  be  those  for  which  the 
asymptotic  decay  is  18  dB/octave  or  better.  This  holds  regardless  of  the  size  of  the 
near-by  sidelobes.  Even  the  window  in  figure  15.  with  a  peak  sidelobe  less  than  -98 
dB,  does  not  perform  as  well  as  the  Hanning  window  which  has  a  -31.5  dB  sidelobe; 
in  fact,  the  error  is  about  70  times  larger.  This  numerical  result  is  consistent  with  the 
discussion  surrounding  figure  19  —  we  have  virtually  eliminated  spillover  across  the 
f  =  0and  f  =  ,5/A  borders  by  using  rapidly  decaying  windows,  when  f(l~.25  A. 

When  the  frequency,  f„,  is  lowered  to  one-tenth  of  the  Nyquist  frequency,  the  best 
windows  are  those  with  decays  of  30  dB/octave  or  better;  see  table  4.  This  trend 
continues  in  table  5  for  f0  equal  to  0.044  of  the  Nyquist  frequency,  where  the  sole 
window  with  a  42  dB/octave  decay  is  observed  to  be  best. 

However,  when  data  frequency,  f,„  approaches  the  bounds  of  (53),  a  different 
ordering  of  windows  is  observed;  see  table  6).  Now  the  windows  with  the  deepest 
near-by  sidelobes  (cases  6  and  15  in  table  2)  have  the  least  error,  although  the  range 
of  errors  in  table  6  are  less  pronounced  than  those  in  the  prev  ious  tables.  Thus,  the 
best  weighting  w(t)  to  employ  depends  on  the  spectrum  of  the  data  sequence  Jx,,} 
that  we  are  subjecting  to  the  DHT.  A  compromise  between  deep  near-by  sidelobes 
and  fast  asymptotic  decay  must  be  accepted  for  a  waveform  occupying  most  of  the 
zero-to-Nyquist  frequency  band.  For  example,  the  average  magnitude-errors  for  the 
minimum  4-term  w indow  of  figure  1 5  are  approximately 

.5E-6,  I.6F-6,  2.6E-6,  6.3E-6  (56) 

for  the  four  frequencies  considered  in  tables  3-6.  This  is  a  very  small  error  and 
indicates  that  the  proposed  DHT  technique  has  considerable  merit,  provided  the 
weighting  in  (50)  is  chosen  appropriately. 

Other  waveforms  were  simulated  and  subjected  to  the  DHT  and  error  analysis. 
They  included  narrow-band  waveforms,  high-frequency  tones  (near  Nyquist),  and 
other  phases  than  0  and  -n/2.  No  surprises  were  encountered;  results  were  con¬ 
sistent  with  those  already  presented  above.  Also  a  simple  Gaussian  weighting  with 
several  different  standard  deviations  was  tried;  the  best  Gaussian  weighting  occured 
when  it  was  set  equal  to  .0053  at  'he  segment  edges  and  .27  at  the  edges  of  the 
retained  points.  Although  the  error  of  the  DHT  could  be  made  to  be  very  good  for 
tone  frequencies  well-removed  from  zero  and  Nyquist  frequency,  it  was  poorer  for 
frequencies  fairly  near  zero  and  Nyquist;  the  worst  error  was  I.36E-4  at  f„A 
0.044,  which  is  20  times  worse  than  can  be  realized  with  weighting  1 5  in  table  6.  This 
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is  due  to  the  spillover  caused  by  the  wide  mainlobe  of  the  Gaussian  window  when 
the  effective  duration  of  the  weighting  is  chosen  narrower  than  the  segment  length. 


Table  3.  Average  Magnitude-Error  for  f0A  =  V^/5  =  .283 


Weighting 

Error  X 

Cosine 

109  for 

Sine 

1 

2.83 

8.52 

2 

1.89 

4.18 

3 

6050. 

7470. 

4 

4350. 

5370. 

5 

10500. 

12900. 

6 

70.4 

86.4 

7 

2090. 

2580. 

8 

1.38 

1.45 

9 

1.75 

3.47 

10 

1.38 

1.45 

1 1 

1.38 

1.45 

12 

1.41 

1.72 

13 

52400. 

64700. 

14 

4700. 

5810. 

15 

415. 

512. 

Table  4.  Average  Magnilude-Krror  for  fyA  =  \^/27  =  .052 


Weighting 

Error  x 

Cosine 

109  for 

Sine 

1 

305. 

454. 

2 

138. 

206. 

3 

16300. 

24400. 

4 

11700. 

17500. 

5 

28200. 

42300. 

6 

164. 

245. 

7 

5580, 

8360. 

8 

.56 

.77 

9 

108. 

161. 

10 

.23 

.23 

11 

.37 

.48 

12 

24.3 

36.2 

13 

142000. 

213000. 

14 

12700. 

19000. 

15 

1090. 

1630. 
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Table  5.  Average  Magnilude-Krror  for  f0A 

=  V^/65  =  .022 

Krror  x  109  for 

Weighting 

Cosine 

Sine 

1 

4360. 

6020. 

2 

1950. 

2700. 

3 

33100. 

45800. 

4 

23300. 

32200. 

5 

58000. 

80300. 

6 

61.1 

52.0 

i 

10700. 

14900. 

8 

39.5 

54.3 

9 

1520. 

2110. 

10 

.93 

1.27 

11 

20.4 

28.1 

12 

326. 

452. 

13 

300000. 

416000. 

14 

25300. 

35100. 

15 

1890. 

2620. 

Table  6.  Average  Magnitude-Krror  for  f0A 

=.  vi/350  =  .004 

Krror  x 

106  for 

Weighting 

Cosine 

Sine 

1 

495. 

501. 

2 

154. 

86.1 

3 

24.0 

50.6 

4 

45.2 

50.5 

5 

61.5 

129. 

6 

7.01 

7.27 

7 

29.3 

15.4 

8 

132. 

241. 

9 

95.4 

45.6 

10 

189. 

398. 

1  1 

34.0 

70.6 

12 

9.07 

15.3 

13 

899. 

429. 

14 

41. 

53.3 

15 

3.7! 

6.31 
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SUMMARY 

When  strong  interference,  either  tonal  or  narrowband,  occurs  additively  with  a 
desired  signal,  its  effect  on  frequencies  removed  from  the  interference  band  can  be 
greatly  reduced  by  using  windows  with  low  sidelobes  and  significant  decay  of  the 
sidelobes.  Thus  close-by  interference  rejection  requires  the  immediate  sidelobe 
region  of  the  window  to  be  small,  while  distant  interference  rejection  requires  a 
rapidly-decaying  sidelobe  response.  The  type  of  windows  considered  here  furnish 
several  alternative  choices,  depending  on  the  application  of  interest,  and  range  from 
-31  dB  to  -98  dB  for  the  peak  sidelobe,  or  6  dB/octave  to  42  dB/octave  for  the 
asymptotic  decay.  The  weighting  given  by  (3)  is  non-negative  for  all  the  numerical 
coefficients  listed  here. 

Use  of  some  of  these  windows  in  a  proposed  discrete  Hilbert  transform  technique 
yields  very  small  errors  over  a  wide  range  of  frequency  components  of  the  given 
data  sequence,  provided  that  50%  overlapped  processing  is  employed.  Two  DFTs 
are  necessary  per  segment,  and  the  edge  50%  of  the  output  values  of  each  segment 
must  be  discarded.  Larger-size  DFTs  yield  less  error,  especially  for  frequencies  near 
the  bounds  in  (53),  but  they  require  more  stc-age. 
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Appendix 

PROGRAM  FOR  DISCRETE  HILBERT  TRANSFORM 

The  subroutine  for  the  DHT  technique  described  earlier  is  contained  in  lines  660- 
820  of  the  program  listed  below  in  table  A~l.  The  remainder  of  the  program  con¬ 
tains  the  data  creation,  the  weightings  in  lines  100-240,  the  error  calculation,  and  a 
DFT  in  the  subroutine  beginning  with  line  830.  On  input  to  the  DHT  subroutine  call 
in  line  500,  array  X  contains  the  N-point  data  segment  to  be  Hilbert  transformed, 
and  array  Y  is  zero.  On  output  from  the  DHT  subroutine  call,  elements  N/4  to 
3N/4  of  arrays  X  and  Y  contain  the  central  half  of  the  original  data  segment  and  its 
DHT,  respectively. 


Table  A-l.  Program  for  Discrete  Hilbert  Transform 


10 

20 

30 

40 

50 

60 

70 

80 

90 

108 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

248 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 


F0de  1  =SQR<2)  -'5  !  5,27,65.350 

Phas e  =  -P I ^2 

N= 1 024  !  SIZE  OF  FFT 

PRINTER  IS  0 

PRINT  "F0de i  »";F0del , "Phase  *";Phase,"H  = “ ; H 
PRINT 

OPTION  BASE  1 

DIM  X<  1024)  ,  V'  1024;- ,  CC  1024),  S<  1024  >  ,  W<  1024) 
DIM  H<0:59>,fl<0:3> 

DATA  .5, .5,0,0 

DATA  .42, .5, .08,0 

DATA  .42659, .49656, .07685,0 

DATA  .42323, .49755, .07922.0 

DATA  . 44959, . 49364, . 05677, 0 

DATA  .35875, .48829, . 14128, .01 168 

DATA  .40217, .49703, .09892, .00188 

DATA  .375, .5, . 125,0 

DATA  .40897, .5, .09103,0 

DATA  . 3125, .46875, . 1875, .03125 

DATA  . 338946, . 431973, . 161054, .018027 

DATA  . 355768, . 487396, . 144232, .012604 

DATA  .53836, .46164,0,0 

DATA  . 4243801 , . 4973406, . 0782793,0 

DATA  . 3635819, . 4891775, .  1365995,  .01064:1 

READ  H<*> 

T  =  2*P I *F0de 1 
FOR  K* 1  TO  N 
B=T*K+Phase 
C < K>  =C0S <  B ) 

SOO*SIN<B) 

NEXT  K 

FOR  Nc  =  l  TO  15 
FOR  L*0  TO  3 
AsL)*H<L+4*<Nc-l >) 

NEXT  L 

PRINT  "<■ .  "  ;Nc  ;  "  >  AC0) ;  A<  1  >;  A<2  «;  At  3> 

A( 1 >»-A< 1 ) 


31 


TR  6239 


380  A<3>=-A<3> 

390  T  =  2*P I / <  N+ 1 > 

400  FOR  K= 1  TO  N 

410  B  =  T  *  K 

420  S  =  Fk0> 

430  FOR  L= 1  TO  3 

440  S=S+A<L>*CGS<L*B> 

450  NEXT  L 

480  W  <  K >  =S 

470  NEXT  K 

480  MAT  X=C 

490  MAT  Y=ZER 

500  CALL  Dht  <  N , X <  *  > , Y <  *  > , W<  *  >  ) 

510  Mr =H i =Ar=A i =0 

520  FOR  k  =  N r 4  TO  3* Nr 4 

530  Dr*ABS<C<K>-Xao> 

540  B  i  =  AES  03  <K  > -Y  <  K  )  ) 

550  Mr  =  MAX  <  Mr  ,  Dr  ') 

580  M i =  M  A X < M i ,  D  1 > 

570  Ar =Ar+Dr 

580  A i =A i +D i 

590  NEXT  K 

80W  PRINT  Mr  ,  M  i  ,  Ar  •'  (.  5  *  N  +  1  )  ,  H  i  /  (  .  5*N  +  1  > 

810  PRINT 

820  NEXT  He 

830  PRINT  LIN'"  3) 

840  PRINTER  IS  18 

850  END 

8  8  W  U  £  D  h  t,  1  N  ,  X  <.  *  >  ,  Y  '■  *  1 ,  M  *  )  >  D I  S  C  R  E  T  E  H  I  L  £  EF'T  T  R  A  N  S  F  0  R  M 

670  MAT  X=X . W 

680  CALL  Fftl0<N,X<*>, Y<*>> 

690  FOR  K  =  2  TO  Nr2 

700  X<K>=X<K)*2 

710  YOO=-YOO*2 

720  NEXT  K 

730  FOR  K=Nr2+2  TO  N 

740  X  <  K>  =  Y  <  K>  =0 

750  NEXT  K 

760  CALL  Fftl0<N,X<*>,Y<*>) 

770  FOR  K=Nr4  TO  3*N/4 

780  T=lr<W<K>*N> 

790  X<K)=XOO*T 

800  Y<K>*-Y(IO*T 

810  NEXT  K 

820  SUBEND 

830  SUB  Fftl0<N,XO),Y(»>>  !  N  <=  2'10  =  1024,  N  =  2  '  INTEGER 
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